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ABSTRACT 

Stationary magnetohydrodynamic (MHD) outflows from a rotating accretion disk are investigated 
numerically by time-dependent axisymmetric simulations. The initial magnetic field is taken to be 
a split-monopole poloidal field configuration frozen into the disk. The disk is treated as a perfectly 
■ conducting, time-independent density boundary [p(r)] in Keplerian rotation. The outflow velocity from 

t^J" ' this surface is not specified but rather is determined self-consistently from the MHD equations. The 

temperature of the matter outflowing from the disk is small in the region where the magnetic field is 
inclined away from the symmetry axis (c 2 <C v k)> but relatively high (c 2 < vj^) at very small radii 
in the disk where the magnetic field is not inclined away from the axis. We have found a large class 
of stationary MHD winds. Within the simulation region, the outflow accelerates from thermal velocity 



OO 
OS 



(~ c s ) to a much larger asymptotic poloidal flow velocity of the order of one- half yX?M/r^ where M is 
the mass of the central object and fj is the inner radius of the disk. This asymptotic velocity is much 
' larger than the local escape speed and is larger than fast magnetosonic speed by a factor of ~ 1.75. 

The acceleration distance for the outflow, over which the flow accelerates from ~ to, say, 90% of the 
asymptotic speed, occurs at a flow distance of about 80ri. The outflows are approximately spherical, 
with only small collimation within the simulation region. The collimation distance over which the flow 
becomes collimated (with divergence less than, say, 10°) is much larger than the size of our simulation 
region. Close to the disk the outflow is driven by the centrifugal force while at all larger distances the 
flow is driven by the magnetic force which is proportional to —W(rB^) 2 , where is the toroidal field. 

Our stationary numerical solutions allow us (1) to compare the results with MHD theory of stationary 
flows, (2) to investigate the influence of different outer boundary conditions on the flows, and (3) to 
investigate the influence of the shape of the simulation region on the flows. Different comparisons were 
made with the theory. The ideal MHD integrals of motion (constants on flux surfaces) were calculated 
along magnetic field lines and were shown to be constants with accuracy 5 — 15%. Other characteristics 
of the numerical solutions were compared with the theory, including conditions at the Alfven surface. 

Different outer boundary conditions on the toroidal component of the magnetic field were investigated. 
We conclude that the commonly used "free" boundary condition on the toroidal field leads to artificial 
magnetic forces on the outer boundaries, which can significantly influence to the calculated flows. New 
outer boundary conditions are proposed and investigated which do not give artifical forces. 

We show that simulated flows may depend on the shape of the simulation region. Namely, if the 
simulation region is elongated in the z— direction, then Mach cones on the outer cylindrical boundary may 
be partially directed inside the simulation region. Because of this, the boundary can have an artificial 
influence on the calculated flow. This effect is reduced if the computational region is approximately 
square or if it is spherical. Simulations of MHD outflows with an elongated computational region can 
lead to artificial collimation of the flow. 

Subject headings: jets, accretion disks — outflows: jets — galaxies: magnetic fields — plasmas — stars 



1. INTRODUCTION 

The existence and nature of magnetohydrodynamic 
(MHD) outflows from an accretion disk threaded by an 
ordered magnetic field is a long-standing astrophysical 
problem. The problem has been investigated theoretically 



by many authors (Blandford & Payne 1982; Pudritz & 
Norman 1986; Sakurai 1987; Koupelis & Van Horn 1989; 
Lovelace, Berk & Contopoulos 1991; Pcllcticr & Pudritz 
1992; Konigl & Ruden 1993; Cao & Spruit 1994; Con- 
topoulos & Lovelace 1994; Contopoulos 1995; Ostriker 
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1997). See also reviews by Bisnovatyi-Kogan (1993) and 
Livio (1997). From the theory, a necessary condition 
for magnetically/centrifugally driven outflows is that the 
poloidal magnetic field at the disk's surface be inclined 
away from the symmetry axis (z) at a sufficiently large 
angle. 

However, the analytical theory makes drastic simplifi- 
cations such as assuming self-similar dependences on the 
radial distance (r in cylindrical coordinates), or by inte- 
grating over the cross-section of the outflow. The self- 
similar solutions have divergences at both small and large 
r so that the influence of these regions is unknown. 

Numerical MHD simulations are essential to establish 
the existence and understand the nature of magneti- 
cally/centrifugally driven outflows. Stationary and non- 
stationary MHD flows were investigated by Kudoh & Shi- 
bata (1995, 1997a, b) in one-dimensional (1.5D) simula- 
tions. These simulations allowed an investigation of out- 
flows for a wide range of parameters. However, they sup- 
posed a fixed configuration of the poloidal magnetic field. 
Two-dimensional (2.5D) simulations of outflows from ac- 
cretion disks were performed by Uchida & Shibata (1985), 
Shibata & Uchida (1986), Stone & Norman (1994), Mat- 
sumoto et al. (1996). These simulations led to strongly 
non-stationary accretion and outflows from the disk. In 
most of these studies, the non-stationarity of the solutions 
is due to the start up conditions with the disk rotating 
but the corona of the disk not rotating. In other cases the 
non-stationarity is due to the disk rotating at a signifi- 
cantly sub-Keplerian rate. These simulations are valuable 
in showing that temporary MHD outflows are possible, but 
the results depend strongly on the assumed initial condi- 
tions. 

In order to avoid the strong dependence on initial con- 
ditions and the problems associated with following the in- 
ternal dynamics of the accretion disk, we earlier proposed 
treating the outer, surface layers of the disk as a boundary 
condition (Ustyugova et al. 1995; Koldoba et al. 1996; 
Romanova et al. 1997; Romanova et al. 1998). This 
approach has been followed by others (Ouyed & Pudritz 
1997; Ouyed, Pudritz & Stone 1997; Meier et al. 1997). 
In these simulations the "disk" represents an outer layer 
of the accretion disk. In actual situation, the outflowing 
matter will affect the disk evolution, or at least to the 
evolution of the surface layers of the disk. The angular 
momentum carried away by MHD outflows can give a disk 
accretion rate much larger than the viscous accretion rate 
of say an a-disk, but the accretion speeds are typically 
much smaller than the free-fall speed (Lovelace, Romanova 
& Newman 1994; Lovelace, Newman & Romanova 1997). 
Thus, the disk can be treated as stationary during the for- 
mation and establishment of MHD outflows which takes 
place on a free-fall time scale. However, the long-time 
simulations of outflows including the back reaction on the 
disk are clearly of interest for future research. 

Different initial magnetic field configurations have been 
assumed in earlier studies. The initial field assumed by 
Ouyed & Pudritz (1997) was the Cao & Spruit (1994) 
field which decreases slowly with radial distance on the 
disk's surface. On the other hand, the initial magnetic field 
of Ustyugova et al. (1995) was the split-monopole field 
(Sakurai 1978; 1985), which decreases rapidly with radial 



distance on the disk surface. The temperature of matter 
outflowing from the disk of Ouyed & Pudritz (1997) was 
small, and the initial magnetic field was weak. However, 
Ouyed & Pudritz (1997) introduced a spectrum of turbu- 
lent Alfven waves with a high pressure which is similar to 
having a high temperature corona. Thus the approach of 
Ouyed & Pudritz (1997) is similar to that of Ustyugova et 
al. (1995) where the magnetic field is weak and the coronal 
temperature is high. In both papers, the initial twist of 
the magnetic field results from the disk rotation because 
the corona is not rotating. This twisting of the magnetic 
field gives the collimation observed in both papers. 

It is important to get stationary outflows using time- 
dependent MHD equations because the non-stationary 
flows may be artifacts of the initial conditions. Station- 
ary magneto-centrifugally driven outflows for relatively 
low temperature of the "disk" matter were obtained in 
the 2.5D simulations by Romanova et al. (1997) for the 
case where the initial magnetic field was a "tapered" split 
monopole type field. This work found that in the station- 
ary state the outflow was quasi-spherical with essentially 
no collimation within the simulation region. Close to the 
disk the outflow was driven by the centrifugal force while 
at larger distances the magnetic force was dominant. 

In this work we investigate the case of a pure (that is, 
non-tapered) split-monopole magnetic field by axisymmet- 
ric (2.5D) numerical simulations. The motivation was to 
study MHD outflows from a relatively cold accretion disk 
where magnetic field lines are inclined away from the sym- 
metry axis. To remove the influence of the region near 
the axis where magnetic field lines are not significantly in- 
clined, we pushed hot matter from the disk in the small 
area around the axis. We compare our simulation results 
with the theory of stationary MHD flows. Further, we 
use our stationary simulation flows to investigate the in- 
fluence of outer boundary conditions. Our earlier study 
(Romanova et al. 1997) showed that some simple outer 
boundary conditions on the toroidal magnetic field can 
lead to artificial forces on the boundary which significantly 
influence the flow within the simulation region. Here, we 
consider in further detail the influence of outer boundary 
conditions on the calculated flows. 

In §2 the theory of stationary MHD flows is briefly re- 
viewed. In §3 the numerical model is presented. The in- 
fluence of the outer boundary condition on the toroidal 
magnetic field and the shape of the computational region 
is analyzed in §4. In §5 we present results of simulations 
of stationary flows and compare them with theory In §6 
conclusions of this work are summarized. 

2. THEORY OF STATIONARY MHD FLOWS 

The theory of stationary, axisymmetric, ideal MHD 
flows was developed by Chandrasckhar (1956), Woltjer 
(1959), Mestel (1961), Kulikovskyi & Lyubimov (1962), 
and others. Under these conditions the MHD equations 
can be reduced to a single equation for the "flux function" 
^(r, z) in cylindrical (r,(j),z) coordinates (Hcincmann & 
Olbcrt 1978; Lovelace et al. 1986). The flux function 4» 
labels flux surfaces so that ^(r, z) =const represents the 
poloidal projection of a field line. The equation for *S> 
is commonly referred to as the Grad-Shafranov equation 
(Lovelace et al. 1986). 
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2.1. Integrals of Motion 

For axisymmetric conditions the flow field can be writ- 
ten as v = v p + w^e^ where v p is the poloidal (r, z) com- 
ponent, «0 = lot > is the toroidal component, and 
is the unit toroidal vector. Similarly, the magnetic field 
can be written as B = B p + B^e^. The ideal MHD equa- 
tions then imply that certain quantities are constants on 
any given flux surface ^(r, z) =const or equivalently they 
are constants along any given stream line or a given mag- 
netic field line. These integrals are functions of \I/ (see for 
example Lovelace et al. 1986), 



Vp 4np Bp 


(1) 


UJr 2_^±^ A(v]/) 


(2) 


w f^=0(*) 
4npr 


(3) 


S = S(V) 


(4) 



w + $ 



1 



fi)V = E(V) (5) 



Here, S is the entropy, w is the enthalpy, and $ is the grav- 
itational potential. The quantity K corresponds to the 
conservation of mass along a streamline, A to the conser- 
vation of angular momentum, £1 to the conservation of he- 
licity, S to the conservation of entropy, and E (Bernoulli's 
constant) to the conservation of energy. 

The remaining MHD equation (which cannot be written 
in the integral form) is the Euler force equation across the 
poloidal magnetic field line (Bogovalov 1997), 



2 89 cos 9 2 2 \ 
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(6) 



which is equivalent to the Grad-Shafranov equation. Here, 
d/dn is the derivative in the direction perpendicular to 
magnetic field lines and directed outward from the axis, 9 
is the angle of inclination of the poloidal magnetic field 
line away from the z— axis, s is the distance from the 
disk along a magnetic field line, va p = |B p |/^/47rp and 
va4> = \B ( f > \/yJ4irp are the poloidal and azimuthal Alfven 
velocities. The quantity 89 /ds is the curvature of mag- 
netic field line. The first two terms in equation (6) are de- 
termined by the non-diagonal (tension) part of the stress 
tensor, pViVk + (p + B 2 /87r) Sik — BiBk/An. The third term 
is determined by the total (matter plus magnetic) pressure 
p + B 2 /87r and the gravity force <9$/<9n. 

2.2. Physical Sense of Integrals of Motion 

To clarify the physical sense of the integrals of motion, 
it is useful to derive the fluxes of mass, angular momen- 
tum (about the z— axis), and energy The corresponding 
conservation laws for stationary conditions are 



V • (pw p ) = , 



(7) 
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B p (vB) 
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0. (9) 



Because v„ 



Bp, the vector flux densities are directed 



along the field lines. 

Consider the fluxes through an annular region with sur- 
face area element dS. The matter flux through the ax- 
isymmetric surface S extending out from the z— axis is 



T M = J dS ■ pw v = ^JdS ■ B P K(V) , 



(10) 



where we took into account the integral (1). dS ■ B p is 
the magnetic flux through the annular region bounded 
by flux surfaces 'J and * + d&. Thus we can change 
from space integration to integration over 'J. Because 
B r = -(l/r)d&/dz and B z = (l/r)dV/dr, we have 



FmW = \j 



(11) 



where ^> — corresponds to the z— axis. Similarly, 



■FlW — J dS ■ (rpVpVj, 



B p ri?0 

4lT 



= \J d9'A(9')K(*') , (12) 



dS ■ pv p 



v 2 B 2 

— + - — + $ 

2 4np 



BpV B 

47T 



= \f d^'{E + ACl)K . (13) 



Thus, Kdfy/2 is the matter flux between the flux surfaces 
separated by d'J/, AKd^/2 is the angular momentum flux, 
and (E + AQ)Kd^/2 is the energy flux. Note that A(f) 
is specific angular momentum carried along the magnetic 
field line = const, E(^>) + Af2(^) is the specific energy, 
and fi(^) is the angular velocity of the disk at the point 
where the magnetic field line or flux surface * = const 
intersects the disk (for |v p | ^ at the disk). 

2.3. Conditions at the Alfven Surface 

Conditions at the Alfven surface are known to be im- 
portant for the global properties of MHD flows (Weber 
& Davis 1967). Equations (2) and (3) constitute a linear 
system of equations for u> and B^. The determinant of 
this system is zero if K 2 = 4irp. Under this condition a 
solution exists if A = r 2 fi (Weber & Davis 1967) which 
corresponds to v p = B p / \J4ixp = va p - This is the con- 
dition which defines the Alfven surface. Figure 1 shows a 
possible field line $ = const and Alfven surface A. The 
radius at which this field line intersects the disk is ^('J'). 
The radius at which it crosses the Alfven surface is 7\4(\P). 
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Fig. 1. — The figure shows a poloidal magnetic field line \P = const (solid line) which starts at the radial distance 
r = rd(^) on the disk and crosses the Alfven surface A (dashed line) at the radial distance ta(^) from the axis. 



The density at this point on the Alfven surface is pa{^)- 
Thus, 



p A (V) = K 2 (*)/4ir , r^(¥)=A(¥)/fi(tt). (14) 
Equations (2) and (3) give 

1 - PAr\l pr 2 



= Q 



1 - Pa/p 



= rQy/AirpA — 



1 - r\/r 2 



Pa/p 



(15) 



(16) 



Taking into account (14) - (16), one can express the fluxes 
of mass, angular momentum and energy, using only the 
values of physical quantities on the Alfven surface: 



^mW = J cM'^Wp-J , 



o 

^(*) = J d^'^pl {E + n 2 r 2 A ) 



(17) 



(18) 



(19) 



2.4. Forces 

For understanding the plasma acceleration, we project 
the different forces onto the poloidal magnetic field lines. 
As mentioned, in a stationary state, matter flows along 
the poloidal magnetic field lines. The acceleration in the 
poloidal (r, z) plane is 



(v p • V)v p + ^(e^ • VX^e,/,) 



(20) 



The last term represents the centrifugal acceleration 
— (u|/r)e r = -rw 2 e r . To get the force per unit mass along 
a magnetic field line, we multiply the Euler equation by a 
unit vector b parallel to B p . This gives 

1 dv d$ 1 - 
f = u; 2 r S m6--^-^ r + — b. VxB xB , (21) 
p os os 4np 

where 9 is the inclination angle of the field line to the 
z— axis. The final term of (21) is the projection of the 



magnetic force in the direction of b, which can be trans- 
formed to 



M 



1 

47T/9 



b- [(V x B) x B] 



1 djrB^f 
Sitpr 2 ds 



which is useful for understanding our results. 

When magnetic field lines are inclined outward, away 
from the symmetry axis, the gravitational force fa = 
d^/ds opposes the outflow of matter from the disk. If 
the matter is relatively cold then the pressure gradient 
force fp = — (1/ p)(dp/ds) is unimportant. Then mat- 
ter can be accelerated outward by the centrifugal force 
fc = uj 2 rsm8 and/or the magnetic force /m- This de- 
termines the driving mechanisms of the outflow, centrifu- 
gal and/or magnetic. The centrifugal force always acts to 
accelerates matter outward if the distance between mag- 
netic field line and the axis increases. Consider the direc- 
tion of the magnetic force. Note that the lines on which 
r£?0 =const are also poloidal current-density lines; that is, 
j p = — (e^/r) x V(rB$) so that j p • V(rB^) = 0. Consider 
a configuration of magnetic field line ^ =const and a line 
of current-density j p as shown in Figure 2. The poloidal 
component of the magnetic force oc — V(r_B^) 2 is perpen- 
dicular to the current-density j p and is shown on Figure 
2 by arrows. Projection of this force onto the poloidal 
magnetic field shows that the force pushes matter upward 
near the disk (lower part of the region), and pushes matter 
downward further from the disk (upper part of the region). 
The (^—component of the magnetic force oc j p x B p acts 
in the direction of the disk rotation and leads to winding 
of the magnetic field line close to the disk and leads to 
unwinding of magnetic field line farther from the disk. 

Thus, magnetic and centrifugal forces may both acceler- 
ate matter, but this depends on the configuration of mag- 
netic field and current-density lines. 

2.5. Collimation 

Consider now the collimation of the flow. From equation 
(6), taking into account that cos# = dr/dn, we have 



(V; 



J Ap 



de 

ds 



1 



Sitpr 2 dn 



8 {rB,f 



cos By 2 




(22) 



At large distances from the Alfven surface r 3> ta, the 
density p < p^, but values pr 2 and v 2 remain finite (Hey- 

vaerts and Norman 1989). Then v 2 Ap = (p/pa) 2 v 2 -c v 2 , 
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Fig. 2. — The sketch shows the directions of the magnetic forces for different configurations of the poloidal magnetic 
field Bp (bold line) and the poloidal current-density j p (thin line). The magnetic force in the poloidal plane oc —W(rB^,) 2 
is perpendicular to the poloidal current-density (bold arrows). The lower part of the figure shows the case where the 
magnetic force pushes matter away from the disk, while the top part of the figure shows the opposite situation. The force 
in the azimuthal direction j p x B^/c acts in the direction of rotation of the disk in the bottom part of the figure line and 
in opposite direction in the top part. 



so that the second term on the left-hand side of (22) is 
negligible. On the right-hand side, only the first term is 
important for r»r4. Then, equation (22) simplifies to 



36 
ds 



1 



d 



8npr 2 dn 



(23) 



Thus, the curvature of magnetic field lines in the region 
r 3> ta is determined by the gradient (rB^) 2 . Figure 3 
shows examples where two magnetic field lines and ^ 
cross the line rB^ = const. Again, the magnetic force 
oc —\I{rB^) 2 acts in the direction perpendicular to the 
current-density j p . Here, we are interested in the projec- 
tion of this force onto a poloidal magnetic field line. From 
the figure one can see that the magnetic force acts to "col- 
limate" the magnetic field line Vf^ and "anticollimate" the 
field line * 2 - 



3. NUMERICAL SIMULATIONS OF MHD OUTFLOWS 

For our time-dependent simulations of axisymmetric 
flows of an ideal plasma in a gravitational field the equa- 
tions are 

'''' ' V • (pv) - , (24) 



at 

d(pv) 
dt 



v • r = Pg 



— - Vx(vxB) = 

dt x ' 



d(pS) 
dt 



V • (pvS*) = 



(25) 
(26) 
(27) 



Here, S is entropy, T jk = pv 3 v k + pSjk + (B 2 <5 jfe /2- 
BjBk) /(47r); is the stress tensor; g = — V<I> is the grav- 
itational acceleration; and $ is the gravitational potential 
of the central object. 

The energy equation (27) is written in conservative 
form. From (24) and (27), one also has 



d[ P f(S)} 
dt 



+ V • [p/(5)v] = , 



(28) 



for any continuous function f(S). We take the equation 
of state to be p — (7 — I )pe, where e is specific internal 
energy, and 7 = const. In the present work 7 = 5/3. We 



take f(S) = pjp 1 ■, because the right-hand side depends 
only on the entropy. We solve the system of equations 
(24)-(26) and (28) numerically. 

The central object of mass M is at the center of our coor- 
dinate system. The disk is located at z = and is treated 
as a perfectly conducting surface rotating with Keplerian 
velocity Vk(t). The gravitational acceleration g diverges 
as r — » 0, but of course the presence of a star or black 
hole changes this dependence. Instead of including the fi- 
nite size of the central object, the gravitational potential is 
smoothed close to the origin, $ = —GM/(r 2 + z 2 + r 2 ) 1 ^ 2 , 
where Ti is the smoothing radius. The value r$ is always 
much smaller than the size of the computational region. 
For this smoothed potential, the Keplerian velocity (for 
z = 0) becomes vk — r\J GM / (r 2 + ) 3 / 4 . Our results 
do not depend significantly on r, because the main part of 
the outflow occurs from the inclined magnetic field in the 
region of the disk where r > r^. 

3.1. Numerical Method 

Equations (24) - (26) and (28) were solved with our Go- 
dunov type numerical code (Koldoba et al. 1992; Koldoba 
& Ustyugova 1994; Ustyugova et al. 1995). The code 
is based on the ideas of Roe (1986) for hydrodynamics 
and the related ideas of Brio & Wu (1988) for MHD. This 
type of TVD numerical scheme has also been developed 
and investigated by others (for example, Ryu, Jones & 
Frank 1995). The code has passed a number of essen- 
tial tests ( Koldoba et al. 1992; Koldoba & Ustyugova 
1994 ) which are analogous to those described by Ryu 
et al. (1995). Compared with our earlier applications of 
this code (Ustyugova et al. 1995; Koldoba et al. 1996), 
a number of improvements have been made, including a 
procedure for guaranteeing that V • B = 0. To satisfy the 
condition V • B = 0, we projected the calculated magnetic 
field to the sub-space of solenoidal functions B at each 
time step. We introduced the function vp, which satisfies 
the equation 



d_ 

dr 



ld_ 

r dr 



dz 2 



* = - r 



dB r dB 2 



dz 



dr 



Then the magnetic field B was calculated for which V-B = 
0. A similar method was used by Ryu et al. (1995). 

Most of the simulations were done on a grid with N r x N z 
points in cylindrical coordinates. For the calculations on 
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Fig. 3. — The figure shows a poloidal current-density line j p on which rB^ = const and two poloidal field lines. On 
the field line the magnetic force oc —V(rB<p) 2 acts to give collimation while for the field line \&2 the force acts to 
"anti-collimate" the flow. 



an approximately square region we used an inhomogeneous 
grid with 100 x 100 points while for the axially elongated 
region we used a homogeneous grid with 50 x 200 points. 
We also did a smaller number of simulations using spheri- 
cal coordinates (R, 9, <j>) and a grid Nr x Ng = 100 x 50. 

3.2. Initial Conditions 

The motivation for this work was the study of station- 
ary MHD flows. Hence it may appear that the initial con- 
ditions are unimportant. However, in practice, an unfa- 
vorable choice of initial conditions can lead to an essen- 
tially longer stage of transition to stationary state, or even 
worse, stationary flows may never be reached. The region 
right above the disk is the most important, because the ve- 
locity distribution near the disk determines the number of 
boundary conditions (the flow may be subsonic or super- 
sonic). Also, the physical parameters, such as density and 
magnetic field, are largest just above the disk. Hence, we 
worked more carefully on the equilibrium at small values of 
z. At large z, approximate equilibrium in the ^-direction 
was sufficient, because the magnetic field, which is dom- 
inant in the corona, stabilizes matter against the violent 
movements. The expressions given below are found to be 
useful initial conditions which give a smooth start up of 
the outflows. 

The initial conditions are arranged as follows. The disk 
and corona are considered to be threaded by a poloidal 
magnetic field of monopole type (Sakurai 1987), B p = 
Q(R- Rq)/|R- Rq| 3 , where Q = B h 2 is the "charge" 
of the monopole, Rq is the position vector of the monopole 
located on the symmetry axis at a distance h below the 
disk. 

The temperature on the disk surface, which is propor- 
tional to the square of (isothermal) sound speed Cy = p/p, 
was taken to have the dependence 



$(r, 0) 



k + K*e 



2 M 



(29) 



where k and are parameters, is a characteristic ra- 
dius inside of which the disk is relatively hot. For speci- 
ficity we take tt = 2rj. For r ^> rx, c^/$(r, 0) w k; 
that is, the sound speed is constant fraction of the Ke- 
plcrian velocity in the region where we expect centrifu- 
gally/magnetically driven outflows from the disk. The 
term in (29) with increases the temperature in the re- 
gion near the axis, where magnetic or centrifugal outflows 
are not expected. For actual conditions, this part of the 
flow may be connected with the star or black hole (Livio 
1997). 



We supposed that the initial temperature of the corona 
is a function only of r, so that the equation (29) is the ini- 
tial condition for the entire computational region. Also, 
we supposed that in z— direction the gravitational force is 
balanced by the pressure gradient, 



ldp 
p dz 



,(r) dp 



p dz 

The solution of this equation is 



dz 



p(r, z) = p d (r) exp 



$(r, 0) - $(r, z) 



(30) 



(31), 



where Pd{f) is the pressure on the disk surface, and $(r, 0) 
is the gravitational potential on the disk. 

In the initial state, the surface of the disk is in equilib- 
rium. We suppose that the gravitational force on the disk 
surface is compensated by the centrifugal force, while the 
matter pressure gradient in r— direction is compensated by 
the magnetic force. That is, on the disk (z = 0), 



dpd 
dr 



1 d{rB^ 
nr 2 dr 



. 



(32) 



Solution of this equation for pressure on the disk pd(f) 
and current Id{r) — (c/2)rB^\ z= Q flowing through a circu- 
lar area of radius r on the disk can be written as 



h{r) 



cBqIi 



Pd(r) = ^cos n 9 



: (l-cos"- 2 <9)-(l-cos" 



-l 1/2 



(33) 



(34) 



where cos(9 = h/(h 2 + r 2 ) 1 / 2 , with h =const, and 9 is the 
inclination of magnetic field line to the axis of rotation. 
This gave us possibility to find initial azimuthal magnetic 
field along the disk, B$(r) = 2Id(r)/cr. Close to the disk, 
rB§ is approximately constant along a magnetic field line. 
The three components of the initial magnetic field on the 
disk are shown in Figure 4. 

To escape rapid twisting of the magnetic field due to the 
difference between the azimuthal velocities of the disk and 
the corona, we supposed that the corona initially rotates 
with an angular velocity which is constant on cylinders 
r =const and equal to vk l r of the disk. As a result of this 
rotation the corona is not in equilibrium in the r-direction. 
However, this lack of initial equilibrium does not disrupt 
the evolution of outflows from the disk, and it does not 
affect the final stationary states where the flow reaches 
equilibrium. 
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Fig. 4.— The figure shows the radial dependences of the different components of the initial magnetic field on the surface 
of the disk. 



3.3. Boundary Conditions 

The lower boundary of our simulation region is the disk 
which is perfectly conducting and rotates at the Kcperian 
rate. Thus the tangential component of the electric field 
in the system of coordinates rotating with the disk is zero, 

[(v - v d ) x B] r , - , (35) 

at z = 0, where = vk&$, and v is the fluid velocity 
just above the disk. This condition means that in this sys- 
tem of coordinates the poloidal velocity is parallel to the 
poloidal magnetic field at z = 0. 

The magnetic field is frozen into the disk so that B z = 
Bdz(r), where B& z is a given function of r and is deter- 
mined by the z-componcnt of the initial monopolc mag- 
netic field. Notice that the two other field components, B r 
and B<p, are not fixed on the disk and change with time 
so as to satisfy the MHD equations in the computational 
region. 

We suppose that the density and entropy on the disk sur- 
face are fixed, p = Pd(f), S = Sd{r) with pd(r) and Sd(r) 
given functions of r which follow from equations (29) and 
(31). 

Note that in the present work, the velocity of outflow 
from the disk is a free variable. This is different from our 
earlier work (Ustyugova et al. 1995). When the velocity of 
outflow from the disk is less than the slow magnetosonic 
speed, then the number of boundary conditions we have 
is sufficient. However, if the outflow velocity is super slow 
magnetosonic, then there should be an additional bound- 
ary condition. Because we do not have this additional 
boundary condition, we suppose that the amplitudes of 
the correponding outgoing waves are equal to zero. This 
is equivalent to the fact that we use the values of calculated 
parameters in the cells just above the disk. 

On the z— axis, all fluxes normal to this axis are equal 
to zero. On the outer boundaries, r = R max or z = Z max , 
the "free" boundary conditions dFj/dn = were used 
for all variables excluding B^. Here, d/dn is the deriva- 
tive perpendicular to the boundary, Fj — {p, f(S), v r , 
V,/,, v z , B r , B z }. Our earlier simulation study (Romanova 
et al. 1997) showed that the condition dB^/dn — can 
lead to unphysical results. The outer boundary condition 
on Bff, is considered in detail in §4. 

4. INFLUENCE OF BOUNDARY CONDITIONS ON FLOW 

If the process of outflow formation is strongly non- 
stationary, then the problem of the influence of outer 
boundary conditions may not appear. This is because it 



is difficult to separate the influence of boundary condi- 
tions from effects connected with non-stationarity. How- 
ever, when the flow goes to a steady-state, we observed 
that the stationary flow pattern can depend on the im- 
posed outer boundary conditions and in some cases on the 
shape of the simulation region. 

It is important to eliminate the influence of boundary 
conditions. It is possible, if (1) the flow is supersonic (su- 
per fast magnetosonic) and it is perpendicular to the outer 
boundaries (then information flows out of the simulation 
region), or (2) the correct boundary conditions are chosen 
by some method. The first condition cannot be realized 
during the stage of establishing of the flow, because ini- 
tially, the flow is subsonic. If the flow is supersonic, but is 
not perpendicular to the boundary, then the Mach cones 
may be partially directed inside the simulation region and 
even supersonic flow may influence the flow inside the re- 
gion. The orientation of the Mach cones depends in general 
on the shape of simulation region. 

The second condition can be realized only in some ap- 
proximation. The "best" outer boundary conditions are 
those which influence only the vicinity of the boundaries 
and not the central part of the simulation region. This 
involves all flow variables, but we will discuss only the 
outer boundary condition on B^,, because we found that 
this condition had the strongest influence on the calculated 
flows. 

The final flow may depend on both the Mach cone ori- 
entation at the boundaries (shape of the region) and on 
the outer boundary condition on B^. In different situa- 
tions one of these factors may be more important than the 
other. To separate their influence on the final flow pat- 
tern, we discuss in §4.1 simulations for a fixed simulation 
region, but with different outer boundary conditions on 
Next, in §4.2 we fixed the boundary condition on B^, 
and investigated the dependence of the flow on the shape 
of simulation region and the Mach cones orientation at the 
outer boundaries. In §4.3 we discuss both factors. 

4.1. Dependence of Flows on Outer Boundary Condition 

on Brj, 

Here, we present results of simulations for a fixed elon- 
gated simulation region R max = 50^, Z max = 200^ 
for three different outer boundary conditions on B^: (1) 
a standard "free" boundary condition, (2) a"force-free" 
boundary condition, and (3) a "force-balance" boundary 
condition. 

4.1.1. 'Free" Boundary Condition 
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First, we performed simulations for the simplest stan- 
dard "free" boundary condition on B^, dB^/dn — 0. We 
observed that this boundary condition may give an ar- 
tificial force on the boundary which influences the flow 
within the computational region. For example, if we sup- 
pose that on the top boundary dB^/dz = 0, then the 
radial component of the current-density equals to zero, 
j r = — (c/A^dB^/dz = 0, which means that the poloidal 
current-density has only a z-component j p = j z z. This 
means that the poloidal current-density j p is not parallel 
to the poloidal magnetic field B p . Consequently, there is 
a force (density) j p x B p /c ^ acting in the <f) direction, 
opposite to the rotation of the disk. Figure 5 shows the 
geometry. 

These 'boundary' forces act such way that the flow never 
reaches a stationary state. To check this fact, and to be 
sure that this is not an effect of non-stationarity of our ini- 
tial configuration, we did simulations for cases which went 
to a stationary state with other outer boundary conditions. 
After establishing stationarity, we substituted the outer 
boundary conditions on B$ to a "free" boundary condi- 
tion. We observed that the stationary state was destroyed 
for the reasons mentioned above. Figures 6a, b demon- 
strate one stage of this destruction, when the poloidal ve- 
locity decreased and became less than fast magnetosonic 
speed in all of the computational region. Even the fluxes of 
mass and other physical parameters through the bound- 
aries are not constants in this simulation. Also, matter 
with magnetic flux enters the region from the right-hand 
side, which is due to the flow being sub-fast magnetosonic. 

To avoid this artificial force, we proposed a "force-free" 
outer boundary condition on B§ (Romanova et al. 1997) 
which we discuss in the next subsection. 

4.1.2. "Force-Free" Boundary Condition 

Another possibility to consider is that the toroidal com- 
ponent of the magnetic force is zero on the outer bound- 
aries. That is, j p || Bp = on the outer boundaries. We 
can write this condition as 

Bp • V(rB ) = . (36) 

We performed simulations with this boundary condition in 
the elongated region and observed that the flow reached a 
stationary state (see Figures 6 c,d). This flow has many 
characteristics of stationary flow. Fluxes of mass, energy, 
and momentum, integrated over different cross-sections, 
are constants. Integrals of motion along magnetic field 
lines are also constants. The flow is well collimated inside 
the simulation region (see Figures 6c, d). However, more 
detailed analysis (see §4.2) shows that this collimation is 
artificial. The "force-free" boundary condition for B^ is 
superior to the "free" boundary condition, because it leads 
to a stationary state, but it does not give the physically 
correct flow. In reality, the magnetic force should not be 
zero on the boundary. There is a magnetic force pushing 
matter outward through the outer boundaries. One can 
see from Figure 6d 

that the poloidal current-density (dashed lines) is not 
parallel to the poloidal magnetic field (solid lines). How- 
ever, on the boundaries (Figure 6d) the two vectors are 
forced to be parallel and thus the poloidal force equals 



zero. This boundary condition is better than the "free" 
boundary condition in the sense that there is no strong 
artificial force at the boundary. From the other side, when 
we put the force equal to zero, it is analogous to applica- 
tion of a force equal to the real force but with the opposite 
sign. 

This is one of the factors which may lead to artificial 
collimation. Another possible factor (Mach cones orienta- 
tion) depends on the shape of the simulation region and is 
discussed in §4.2. 

4.1.3. "Force- Balance" Boundary Condition 

As a next step for improving the outer boundary condi- 
tion on i?^, we take into account the fact that the magnetic 
field is not force-free and j p is not parallel to B p . We start 
from equation (16) for B^ and write it in the form 

rB^ = n^/I^ - f « -SI J— pr 2 , (37) 

1 - pa/p V PA 

where we assume that the density at the boundary is much 
less than that at the Alfven surface, p <C pa for r 2 3> r\. 
Then, we obtain 

Bp • V(rB ) = -n J— Bp • V(pr 2 ) 

V PA 

= rB^Bp ■ W)\n(pr 2 ) = aB r B^ , (38) 

where we supposed that pr 2 = F(^)r a and took into ac- 
count that f2 and pa are constants along magnetic field 
lines. 

Finally, we obtain the outer boundary condition as 

Bp • V(rB ) = aBrBj,, (39) 

where a is a parameter. In this case we got stationary 
flows which are not collimated in the simulation region 
(see Figures 6e, f). Fluxes through the outer surfaces and 
integrals along magnetic field lines are well conserved, as 
in the case of collimated flow, described in §4.1.2. 

The question arises, which boundary condition is cor- 
rect, "force-free" or "force-balance" ? The "force-balance" 
condition is clearly the physical condition because it does 
not generate an artificial force on the boundary. How- 
ever, it is more difficult to apply because there is no di- 
rect method for determining the parameter a. It can only 
be obtained iteratively using additional simulations, which 
is very time consuming. Our analysis indicates that the 
"force-free" boundary condition gives good results as com- 
pared with those obtained using the "force-balance" con- 
dition if the shape of the simulation region is not elongated 
in the z— direction. 

Below, we investigate different runs for "force-free" 
outer boundary conditions on B,p, but for different shapes 
of the simulation region. 

4.2. Dependence of Flows on Shape of Region: 
Orientation of Mach Cones 

We noticed empirically that results of simulations de- 
pend significantly on the shape of simulation region. The 
ratio between R max to Z max is critical. We observed that 
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Fig. 5. — The figure demonstrates the artificial force which can appear on the top outer boundary of the simulation region 
in the case of a "free" boundary condition on B^. Here, B p and j p are the poloidal magnetic field and current-density 
(filled vectors). The hollow arrow shows the artificial force acting in the azimuthal direction. 
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when the region is elongated in z— direction, then the flow 
has tendency to collimate. When the region is square, or 
spherical, or elongated in r— direction, then the outflow 
is almost spherical, that is, only slightly collimatcd. Here, 
we present results of simulations all with "force-free" outer 
boundary conditions but different shapes of the simulation 
region. 

First, we investigated the case where the height of the 
region is the same as before, Z max — 200rj, but the region 
is much wider, R ma x = 170rj. Figure 7 shows that in this 
case we got almost spherical outflow, which is very differ- 
ent from the well-collimated outflow in the narrow region 
at the same boundary conditions (Figures 6c, d). 

We also performed similar simulations in spherical co- 
ordinates with Rmax = 170ft, and got similar result. The 
question is why the flows are so different for different 
shapes of the simulation region? In all cases the flow is 
super fast magnetosonic in most of the region. However, 
note that even if the flow is super fast magnetosonic, infor- 
mation can flow in from the boundaries to the simulation 
region, if the Mach cones are directed inside the simulation 
region. 

The Mach cone projected onto the poloidal plane has a 
half opening angle ip which is 



tan 2 if 



K + c 2 M-v 2 c ) 



-.2 



(40) 



where c sm and c/ m are the slow and fast magnetosonic ve- 
locities, respectively, which satisfy c A s j m — c 2 f m {c 2 s + v 2 A ) + 
= (with v\ = B 2 /(4^p) and v 2 Ap = B 2 /(4^o)) and 



v c = va p c s /{v 2 a + c 2 ) 1 / 2 is the "cusp" velocity (Polovin & 
Demutskii 1980; Lovelace et al. 1986; Bogovalov 1997). 
Figures 6-8 show the Mach cones on the outer boundaries 
for different shapes of the simulation region. We find that 
in the case of an elongated region (Figure 6d) an essen- 
tial part of the Mach cones is directed into the simulation 
region, whereas in the case of an almost square region (Fig- 
ure 7b) only very small part of the Mach cones is directed 
into the region. Figure 8 shows that the most desirable ge- 
ometry - where information flows outward across the outer 
boundary - is obtained in spherical coordinates where all 
Mach cones are directed outward from the simulation re- 
gion. 

The elongated region is the least desirable and as dis- 
cussed it gives artificial collimation of the flow. Note that 
the first stationary MHD flow solutions (Romanova et al. 
1997) were obtained using a simulation region elongated 
in the r-direction. 



4.3. Discussion of Boundaries 

Regarding the outer boundaries, we conclude that sim- 
ulated flows may depend on both the outer boundary con- 
dition on B^ and on the shape of simulation region (the 
Mach cone orientation on the outer boundary) . The influ- 
ence of each of these factors may be different in different 
situations. 

The orientation of the Mach cones at the boundary is 
not connected directly with existence and configuration of 
a stationary flow. However, if Mach cones are partly di- 
rected inside the simulation region along part of the outer 
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Fig. 7. — Results of simulations for an approximately square simulation region R max = 170rj, Z max — 200ri using the 
"force- free" outer boundary condition on B^. Here, r and z are measured in units of rj, which is the inner radius of the 
disk. The lines and arrows have the same meaning as in Figure 6. 




Fig. 8. — Results of simulations for the same case as Figure 7 but using a spherical coordinate system and simulation 
region. 



boundary, then the question arises: what is the result 
of this influence, and how strong is it? Our simulations 
shown that for a "force-free" boundary condition on 
the result is artificial collimation of the flow, whereas in 
the case of a "free" boundary condition there is destruc- 
tion of a stationary flow which was arranged as an initial 
condition. 

From comparison of cases shown in Figures 6d and 7b 
it is not clear that Mach cones are responsible for the col- 
limation of the flow in the case shown in Figure 6d. In the 
narrow region, the magnetic field at the right-hand, outer 
boundary is much stronger than in the case of wide region, 
so that influence of the non-exact "force-free" boundary 
condition should be stronger in the case of a narrow region. 
To check this possibility, we performed simulations in a 
small square region with R max = 50ri and Z max = 50ri 
and found uncollimated almost spherical outflow. This in- 
dicates that the shape of the region is the most important 
factor affecting collimation. 

Another question are evident. Why in the case of the 
"force-balance" boundary condition on B^ in the elon- 
gated region do we find the physically correct stationary 
solution? We suggest that during establishment of the 
stationary flow, which may be quite violent (in spite of 
almost stationary initial conditions), this boundary con- 
ditions kept approach to stationarity less violent (than in 



the case of "force-free" boundary conditions) and kept the 
Mach cones directed outward most of the time. The fact 
that in this case a small part of the Mach cones is directed 
inside the region, means that some inflow of information 
from the outer boundary may have only a small affect on 
the flow. 

For some purposes, such as the study of the propagation 
of jets, it is attractive to use a long narrow computational 
region. The general conclusion of this section is that a nar- 
row region can lead to artificial collimation of the flow or 
invalid solutions unless special care is given to the bound- 
ary condition on B$. 

5. STATIONARY FLOWS: COMPARISONS OF SIMULATIONS 
WITH THEORY 

Here, we describe results of simulations for a region 
Rmax = 170ri and Z max = 200^ with a "force- free" outer 
boundary condition. Simulations of the flow in the same 
region but with the "modified" boundary condition gave 
similar results. 

Figure 9 shows the initial distribution on the disk of the 
Keplerian velocity vk, the poloidal Alfven velocity va p , 
and the sound speed c s . Matter outflowing from the disk 
has a time-independent distribution of density as a func- 
tion of radius. The velocity of outflow from the disk is 
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determined by the solution of the MHD equations in the 
simulation region. The simulations show that in a station- 
ary state matter just above the disk has a velocity some- 
what larger than the slow magnetosonic velocity. This is 
in accord with the theory which indicates that the slow 
magnetosonic surface is located inside the disk (Lovelace, 
Romanova, & Contopoulos 1993). 

Moving away from the disk, matter starts from a low 
velocity, is gradually accelerated and crosses the Alfvcn 
and fast magnetosonic surfaces (see Figure 7a). These 
surfaces are almost parallel to the disk. Figure 10a shows 
the variation of different velocities along a representative 
magnetic field line, the third line away from the z— axis 
in Figure 7, on the flow distance s from the disk. This 
field line, which we refer to as the "reference" field line, 
crosses the top boundary at about the midpoint of this 
boundary. This line is not special, but it is inclined suf- 
ficiently to the axis that magnetic/centrifugal forces are 
important. Figure 10a shows in particular the dependence 
of the poloidal velocity v p (s), which becomes larger than 
the Alfvcn velocity va p fairly close to the disk, at s > lOr^. 
Further, v p becomes larger than the local escape velocity 
v esc for s > 25r^. At larger distances, v p becomes larger 
than the fast magnetosonic velocity c/ m at s > 40ri, and it 
approaches an asymptotic speed which is about 1.75 times 
Cf m at the outer boundary of the simulation region. The 
poloidal velocity is parallel to the poloidal magnetic field 
to a good approximation in accord with the theory. Figure 
10b shows the dependences of v p (s) for different field lines. 

Within the simulation region, the outflow accelerates 
from thermal velocity to a much larger asymptotic poloidal 
flow velocity of the order of Q.hy/GMjrl. Thus, the ac- 
celeration distance for the outflow, over which the flow 
accelerates from ~ to, say, 90% of the asymptotic speed, 
occurs at a flow distance of about 80rj. 

5.1. Mechanism of Acceleration 

Figure 11 shows the different forces acting along the 
"reference" magnetic field line. The centrifugal force (Fq) 
is larger than the magnetic (Fm) or pressure gradient force 
(Fp) immediately above the disk s < 10rj. The magnetic 
force is few times larger than the centrifugal force for larger 
distances, s > 10rj. Note that the pressure gradient force 
is negligibly small. Thus, the main driving forces pushing 
matter outward are magnetic and centrifugal. 

Each poloidal magnetic field line is labeled by its "J/ 
value, which equals the magnetic flux through the circular 
region between the axis and the field line. increases from 
zero on the axis to a largest value on the field line most 
distant from the axis. We integrated the forces to obtain 
the total work performed by the magnetic, centrifugal and 
other forces, along different field lines from the disk to the 
outer boundary. Figure 12 shows the dependence of this 
work on One can see that near the axis (small \P) the 
main work is performed by the centrifugal force, while the 
magnetic force is also important. The work along the "ref- 
erence" field line marked by "i?" on the axis, is done 
mainly by the magnetic force with the centrifugal force 
also important. Going away from the axis to larger \& and 
more inclined magnetic field lines, one can see that the 
magnetic force is more and more important role. Note 
that the work done by the pressure gradient is small for 



all field lines. 

5.2. Analysis of Stationarity 

A first indication of stationarity is when the fluxes of 
mass and other physical quantities become constants in 
time. We observed that the fluxes of mass Tm , angu- 
lar momentum Tl, and energy Te calculated through the 
middle of the region z = 0.5Z max become constants after 
about t > 200ti, where U = 2-KTi/vi where Vi = J GMjri 
and Vi is the inner radius of the disk. The time depen- 
dence of the fluxes is shown in Figure 13. We performed 
simulations for much longer times, t ~ 3500ii, and ob- 
served that these fluxes were accurately constants in time. 
Note that this time corresponds to only t = l.Qt ou t, where 
tout = U (r out /ri) 3 ^ 2 = 2216i,. This indication of sta- 
tionarity is necessary, but not a sufficient sign of a valid 
stationary MHD flow. 

Another indication of stationarity is that the poloi-dal 
velocity becomes parallel to the poloidal magnetic field. 
We observed, that the two vector fields become parallel 
to a high accuracy only after t > t out . Figure 7a shows 
that the two vector fields are close to being parallel even 
at earlier times. 

One can get more complete information about station- 
arity and validity of the MHD flow by comparing the the- 
ory reviewed in §3 with the simulation data. First, the 
integrals of the motion, A, K, E, f2, and S (equations 
4-8) should be constants along any magnetic field line. 
We checked this by numerically calculating these integrals 
along the "reference" magnetic field line. The calculated 
integrals are constants with good accuracy. For exam- 
ple, |Afi|/fi £ 0.06 and \AS\/S <, 0.15. Figure 14 shows 
variation of the integrals. Note, that the integrals are not 
strictly constants in the region immediately above the disk, 
because the grid is not fine enough in this region due to the 
strong gravitational force. Note that the integrals become 
constants as a function of \I> much later (t ^ t out ) than 
fluxes of mass, angular momentum, and energy become 
constants in time. 

Other comparisons of simulations with theory have been 
done. For example, from the theory of stationary flow it 
follows that fluxes of matter, angular momentum, and en- 
ergy flowing inside a given flux tube should be equal to 
fluxes integrated over the Alfvcn surface, equations (17)- 
(19). We calculated these integrals in two ways, using the 
data from our simulations. Figures 15 a, b show these in- 
tegrals as a function of They almost coincide in most of 
the region, excluding the region of large values of \P. The 
latter field lines have a high angle of inclination relative 
to the axis and do not pass through the fast magnetosonic 
surface. These lines are marked by the long-dashed lines 
on one of the curves, and by letter "/" on the axis. 

Figure 16 shows the dependence of the ratio of the 
radii where a magnetic field line crosses the Alfven surface 
and the disk, A = r^*)/^^). This ratio is of interest, 
because the value of the angular momentum per unit mass 
carried by the outflowing matter can be calculated as 

k = Slr 2 A = \ 2 [GMr d {*)\ 1 ' 2 , (42) 

from equation (14). The fact that A is almost constant 
means that the specific angular momentum is proportional 
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Fig. 9.— Dependencies of velocities on r (in units of ri) along the disk in the stationary state, where ri is the inner radius 
of the disk. Here, vk is the Keplerian velocity, c s the sound speed, and va p the poloidal Alfvcn velocity, measured in 
units of Uj. The point 9 = 30° shows the location on the disk where the poloidal magnetic field is inclined to the z— axis 
at an angle 9 = 30°; for larger r we have 9 > 30°. We show only r < 35ri, because the magnetic field lines which start at 
larger r are highly inclined away from the z— axis and also do not pass through the fast magnetosonic surface within our 
simulation region (see Figure 7). Inside the radius tt the disk is relatively hot (see equation 29). 




Fig. 10. — The top panel (a) shows the dependences of different velocities on distance s measured in units of r, from the 
disk along the "reference" magnetic field line. The velocities are measured in units of Vi. The "reference" field line is the 
third field line away from the z— axis in Figure 7a. This field line crosses the top boundary about in the middle. This field 
line "starts" from the disk at r w 6n where it has an angle 9 w 28° relative to the z— axis. Here, v p is the poloidal velocity 
along the field line and v± is the poloidal velocity perpendicular to the field line. Also, va p is the poloidal Alfven velocity, 
Cf m is the fast magnetosonic velocity, and v esc is the local escape velocity. The bottom panel (b) shows the dependences 
of v p (s) for different field lines indentified by the number of the field line counted away from the z - axis in Figure 7a. 




/ 

/ 
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Fig. 11.— Forces acting along the "reference" field line. M is the magnetic force, C the centrifugal force, P the pressure 
gradient force, and G the gravitational force. The distance s is measured in units of n. The scale for forces is arbitrary. 



to y/¥d and can be estimated in this way. Specifically, 
AA/A <J 0.13 for field lines which cross the fast magne- 
tosonic surface. 

The fluxes of mass, energy, and angular momentum flow- 
ing out from the disk depend of course on the magnetic 
field strength on the disk. Figure 17 shows the dependence 
of the matter outflow rate on the disk magnetic field. This 
dependence is analogous to that derived by Kudoh & Shi- 
bata (see Kudoh & Shibata 1995, figure 2b, and Kudoh 
& Shibata 1997b, Figure 24b) who performed 1.5D anal- 



ysis of stationary MHD flows at the fixed configuration of 
poloidal magnetic field. 

5.3. Collimation 

The stationary MHD outflows we find are approximately 
spherical outflows with relatively small collimation within 
the simulation region. Thus, the collimation distance over 
which the flow becomes collimated (with divergence less 
than, say, 10°) is much larger than the size of our sim- 
ulation region. Figure 18 shows the dependence of the 
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Fig. 12. — Work done by the different forces along different field lines from the disk to the outer boundary. Each field 
line is labeled by its value of The field line corresponding to our "reference" line is marked "r," and the "diagonal" 
field line which goes through the top right corner of the simulation region is marked "d." The point "/" is the largest 
radius at which the flow goes through the fast magnetosonic surface. The letters M, C, P, and G stand for magnetic, 
centrifugal, pressure gradient, and gravity. 
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Fig. 13.— Fluxes of mass Tm, angular momentum Tl-, and energy Te across the area z — 0.5Z max as a function of time. 
Time is measured in units of t out ~ 2200ii, where U = 2nri/vi, where r 4 is the inner radius of the disk. 




Fig. 14.— Numerically calculated "integrals of the motion" (equations 1-5) as a function of distance s along the "refer- 
ence" magnetic field line. In ideal MHD, the integrals should be strictly constant. In this plot the scale is such that the 
maximum of the E integral is unity. The other curves have been displaced downward for clarity. 



angle between the poloidal field direction and the z— axis 
on the distance along the "reference" magnetic field line 
(#i ) and along the "diagonal" field line which goes through 
the top right corner of the simulation region (62)- Both 
angles are relatively large (> 30°) near the disk and then 
gradually decrease at larger distances s along the field line. 
This means that some collimation occurs near the disk but 
decreases at larger distances. The angle 6\ for the "refer- 
ence" field line changes from 35° at the disk to 18° at the 
top boundary where the angle of the position vector from 
the origin to the z— axis is about 20.5°. The angle 6*2 for 
the "diagonal" field changes from 49° at the disk to 31° at 
the top boundary where the angle of the position vector 
from the origin to the z— axis is about 39°. 

An important question is whether the outflow becomes 
collimated at large distances to form a cylindrical jet par- 
allel to the rotation axis, or it continues as a wind without 
collimation. To obtain information on this question, we 
calculated the derivatives 89/ ds along magnetic field lines 
as also shown in Figure 18. For both the "reference" and 
the diagonal field lines the derivatives decrease and be- 
come very small at the outer boundary. It appears that the 
derivatives continue to decrease, which would mean that 



the collimation decreases and goes to zero. (The turns at 
the ends of the lines are connected with boundaries and 
do not represent a real collimation effect.) However, the 
present study does not rule out the possible magnetic col- 
limation of the flow at much larger distances. Separate 
simulations in a much larger region are needed to answer 
this important question. 

Earlier, Sakurai (1987) obtained stationary flow solu- 
tions for a split-monopole magnetic field and found flows 
with very gradual collimation at large distances from the 
origin. Our results are similar in this respect to those of 
Sakurai. However, it is not clear that the flows will mag- 
netically collimate to cylinders as predicted by Heyvaerts 
& Norman (1989). Simulations on a much larger region 
arc needed to answer this question. 

It is important to note that analytic, self-similar solu- 
tions for outflows for cases of very gradually decreasing 
poloidal magnetic field in the disk (unlike the present split- 
monopole field) show magnetic collimation with increasing 
distance z from the origin (Contopoulos & Lovelace 1994, 
Contopoulos 1995; Ostriker 1997). Further, note that out- 
flows may be collimated hydrodynamically by the pressure 
of surrounding, ambient matter (Lovelace, Berk, & Con- 
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Fig. 15. — Panel (a) shows the matter flux and panel (b) the angular momentum flux as a function of "J/ calculated in 
two ways. The solid line shows the integrals calculated using equations (11) and (12). The dashed line shows the integrals 
calculated on the Alfven surface using equations (17) and (18). Points "r," "c£," and "/" on the horizontal axis are the 
same as in Figure 12. The long-dashed line shows the region of strongly inclined field lines which do not cross the fast 
magnetosonic surface. These field lines are separated from the other field lines by "/" on the \& axis. 




Fig. 16.— The ratio A = rA(*)/r<j(*) as a function of ^. Points "r," "d," and "/" are the same as in Figures 12. 



v„/v. 



Fig. 17. — The dependence of matter flux from the disk on the ratio va p /vk evaluated at the point r = t\. Here, 
we changed only the magnitude of the magnetic field B p while all other parameters were kept the same. The triangle 
indicates the main case of the paper while the circles were obtained from different runs. The dashed line is simply a curve 
through the points. 



topoulos 1991; Frank & Mellema 1996, Mcllcma & Frank 
1998). This mechanism of collimation needs a separate 
numerical investigation. 

It is of interest to compare the present results on col- 
limation with those of our earlier studies, Ustyugova et 
al. (1995) and Romanova et al. (1997). Ustyugova et al. 
(1995) introduced the treatment of the disk as a boundary 
condition and found non-stationary but well-collimated 
outflows. Ustyugova et al. considered outflows from a hot 
accretion disk where the sound speed c s <J vk and a weak 
magnetic field v\ « v\ on the disk surface. For such 
conditions, the main force driving the outflow is the mat- 
ter pressure gradient while the magnetic force is smaller. 
Also, Ustyugova et al. (1995) used non-equilibrium ini- 
tial conditions where the rotation of the disk was started 
at t = with the corona of the disk not rotating. These 
conditions led to the formation of a strong toroidal mag- 
netic field (by the winding up of the poloidal field) and 



a strong outward propagating torsional Alfven wave. The 
fact, that Alfven velocity was much smaller than Keplerian 
veloicity allowed the build up of the toroidal field which in 
turn gave strong collimation of the outflow. At later times 
in the simulation, the twist of the field relaxed, but the 
matter pressure force continued to push matter along the 
collimated magnetic field lines. Thus, the Ustyugova et al. 
(1995) flows are essentially different from the stationary 
flows discussed in this paper where the dominate driving 
forces are magnetic and centrifugal with the matter pres- 
sure force negligible. 

Romanova et al. (1997) was the first simulation study 
to obtain stationary magneto-centrifugally driven outflows 
with relatively small matter pressure force. The initial 
poloidal magnetic field was a "tapered" monopole config- 
uration. The outflows were found to be uncollimated and 
are therefore similiar to those discussed in this paper for 
the split monopole initial field. 
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Fig. 18.— The figure shows the dependence of the angle between the poloidal magnetic field and the z— axis as a 
function of distance s along the field line. The dependence of the derivatives d6/ds is also shown. The label "1" indicates 
our "reference" magnetic field line and "2" the "diagonal" field line discussed in the text. 



5.4. Illustrative Physical Values 

Here, we discuss physical values of parameters for the 
case of MHD outflows from the disk around a young star. 
The mass of the star is considered to be M = M Q w 
2 x 10 33 gm, and the inner radius of the disk is rj = 10 11 cm, 
which may be somewhat larger than the star's radius. The 
magnetic field threading the accretion disk may arise from 
the "shearing off" and opening of the intrinsic stellar field 
as discussed by Lovelace, Romanova, & Bisnovatyi-Kogan 
(1995). The reference field strength at the inner edge of the 
disk is taken to be Bi = |B p (r i; 0)| = 300 G. The compu- 
tational region extends from (r, z) = to (R m ax, Z max ) = 
(170^,200^) so that R max = 1.7 x 10 13 cm. The char- 
acteristic velocity at the inner edge of the disk is Vi = 
{GM/ri) 1 / 2 w 365 km/s. With our smoothed potential 
the Kcplerian velocity of the disk at r, L is w 0.6uj. 

The characteristic time at the inner edge of the disk is 
ti = 2-KTijvi k> 4.8 hr. The density on the surface of the 
disk at r — ri can be written in terms of the Alfven veloc- 
ity as Pl = B?/(4ttv 2 Ap ) m 6.8 x HT^/tUp) 2 g/cm 3 . 

We find that the mass outflow rate from the top side of 
the disk has the dependence 

M+ = T M r 2 B 2 /v Ki , 

— *Mt^) s (^)W- 

(43) 

where Tm — ^M{vA/vi) is the stationary state value of 
J~m shown in Figure 13 for our reference case. Figure 17 
shows the dependence of Tm on VA P /vi. For the reference 
case, T M ~ 5.6 so that M + w 2.2 x 10~ 6 M /yr. The 
total mass outflow from the disk is evidently 2M+. 

Similarly, we can write the energy outflow rate from the 
top of the disk as 

E+ = T E r 2 B 2 v Kl , 

~ 1 5 x io 34 ^ P (^^\ 3/2 (_*_y (E\ k 

-1.5X10 s ^^ 10 n cm J ^ 300G J [ Mq ) • 

For given distribution of velocities along the disk, and 
given physical parameters at the inner edge of the disk, 
we get T E ~ 1-2 and E + w 2.5 x 10 35 erg/s. The total 
energy outflow from the disk is 2E+. 

Matter is accelerated from v p 0.05ui s=s 18 km/s near 
the surface of the disk to v p w 0.47«i = 170 km/s at the 
outer boundary of the simulation region (see Figure 10). 



6. CONCLUSIONS 

We have studied MHD outflows from a rotating, con- 
ducting accretion disk using axisymmetric simulations The 
disk was treated as a boundary condition, and the initial 
poloidal magnetic field was taken to be a split-monopole. 
The main conclusions of this work are: 

1. In many different runs we observed the formation 
of stationary MHD outflows from the disk. Close 
to the disk the main driving force is the centrifugal 
force. At larger distances the main driving force is 
the magnetic force oc — \7(rB^) 2 . The pressure gra- 
dient force is much smaller than these forces and it 
has no significant role in driving the outflows. 

2. For the considered conditions, the slow magnctosonic 
surface lies inside the disk. Above the disk, the flow 
accelerates and passes through the Alfven and fast 
magnetosonic surfaces, which are almost parallel to 
the disk. Within the simulation region, the outflow 
accelerates from thermal velocity (~ c s ) to a much 
larger asymptotic poloidal flow velocity of the order 
of 0.5^/GM/ri, where M is the mass of the cen- 
tral object, and is the inner radius of the disk. 
This asymptotic velocity is much larger than the lo- 
cal escape speed and is larger than fast magnetosonic 
speed by a factor of — 1.75. The acceleration dis- 
tance for the outflow, over which the flow accelerates 
from — to say 90% of the asymptotic speed, occurs 
at a flow distance — 80r,. 

3. The outflow is only slightly collimatcd within the 
simulation region. The collimation distance for the 
outflow, over which the flow becomes collimated 
(with divergence less than say 10°), is much larger 
than the size of our simulation region. This "poor" 
collimation is similar to that found in our earlier 
work (Romanova et al. 1997) using a different initial 
magnetic field and is qualitatively similar to the very 
gradual collimation found by Sakurai (1987). MHD 
simulations using much larger computatinal regions 
are needed to determine the collimation of the out- 
flow at large distances. Further, separate simulations 
are also needed to study collimating influence of an 
external medium (Lovelace et al. 1991, Mellema & 
Frank 1998). 

4. The stationarity of the MHD flows was checked in a 
number of ways, including comparisons of simulation 
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results with predictions of theory of stationary ax- 
isymmetric flows. We found that: (a) Fluxes of mass, 
angular momentum, and energy across the surface 
z = 0.5Z max become independent of time with high 
precision at early times of simulations t < 0.1t out , 
where t out « 2200i; and t, L = 2-KrJ y/GM/n. (b) 
Integrals of the motion become constants on flux sur- 
faces with accuracy 5% — 15% for t ^ t out . (c) Vec- 
tors of poloidal velocity are parallel to those of the 
poloidal magnetic field lines to a high accuracy. 

5. Different outer boundary conditions on the toroidal 
magnetic field were investigated. We analyzed 
simulation results using and found that collimation 
of the jet and other characteristics of the flow de- 
pend critically on the outer boundary condition on 
i?0 (as well as the shape of the simulation region as 
discussed below) . We observed that the outer "free" 
boundary condition on leads to an artificial force 
which can give apparent magnetic collimation of the 
flow. "Force-free" and "force-balance" outer bound- 
ary conditions were also investigated. The "force- 
free" outer boundary condition was found to give 
valid flow solutions if the simulation region is not 



narrow in r— direction (compared with z-direction). 

6. The question of the optimum shape of simulation re- 
gion was investigated. We have shown that if region 
is narrow in the r— direction, then an essential part 
of the Mach cones on the outer boundaries may be 
directed towards the inside of the computational re- 
gion. This can lead to the influence of the boundary 
on the calculated flow and to artificial collimation. 
This effect is reduced or absent if the computational 
region is approximately square, if it is elongated in 
the r— direction, or if it is spherical. In these cases 
the Mach cones tend to point outside of the com- 
putaional region. 
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